    simpson(y, x, dx, dim, even)

Integrate y(x) using samples along the given axis and the composite
Simpson's rule. If x is None, spacing of dx is assumed.
If there are an even number of samples, N, then there are an odd
number of intervals (N-1), but Simpson's rule requires an even number
of intervals. The parameter 'even' controls how this is handled.

Parameters
----------
y : array
    Array to be integrated.
x : array, optional
    If given, the points at which `y` is sampled.
dx : float, optional
    Spacing of integration points along axis of `x`. Only used when `x` is None. Default is 1.
dim : int, optional
    Axis along which to integrate. Default is the last axis.
even : String["avg", "first", "last"], optional
    'avg' : Average two results: 1) use the first N-2 intervals with
              a trapezoidal rule on the last interval and 2) use the last
              N-2 intervals with a trapezoidal rule on the first interval.
    'first' : Use Simpson's rule for the first N-2 intervals with
            a trapezoidal rule on the last interval.
    'last' : Use Simpson's rule for the last N-2 intervals with a
           trapezoidal rule on the first interval.

Notes
-----
For an odd number of samples that are equally spaced the result is
exact if the function is a polynomial of order 3 or less. If
the samples are not equally spaced, then the result is exact only
if the function is a polynomial of order 2 or less.

Examples
--------

```
>>> x = 0:9
>>> y = 0:9
>>> simpson(x, y)
40.5

>>> y = x .^ 3
>>> simpson(y, x)
1642.5

>>> simpson(y, x, even="first")
1644.5

>>> simpson(y, x, even="last")
1640.5
```